%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------------  Estimate exporter parameters  --------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all

clear 

machine = '';

set(0,'defaultTextInterpreter','latex')

M = 15;

g_euler = 0.57721;

scale = 10^5;

%% Load cost estimates

path_ships = '';
load(path_ships);

%% Compute value for exporters

path_value = '/Users/giulia/Dropbox/giulia-myrto-theo/Data/Entry/Total_Week_Production_in1000MT.csv'; %path to a file containing the average weekly production of raw materials by region.

r_tmp = csvread(path_value,1,1);

r_tmp = r_tmp(:,[1,2,5]);

r_tmp(r_tmp==-100) = nan;

r = zeros(15,15);

expt = r_tmp(:,1);
imp = r_tmp(:,2);

ut_rate = 0.65;

for i = 1:15
    for j = 1:15
        
        if sum((expt==i) & (imp==j))==1 && i~=j
            r(i,j) = r_tmp((expt==i) & (imp==j),3).*dwt(i,j);
        else
            r(i,j) = nan;
        end
    end
end

r = r.*ut_rate.*1000/scale;

%% Exporters' costs
%
% We estimate the exporters' inventory costs and bargaining parameters
% matching simulated and observed trading prices
%
options = optimset('Algorithm','interior-point','display','off','MaxIter',90000,'MaxFunEvals',90000,'TolFun',10^-7,'TolX',10^-7,'TolCon',10^-7);

continent = [1,1,2,2,2,3,4,4,3,5,5,6,7,6,7];
countries = 1:15;

W = Var;
x_exp = zeros(M,2);
meanprices_sim = zeros(M,M);
for i = 1:ng
    n = sum(continent==i);
    ct = countries(continent==i);
    gamma0 = rand(n,1);
    cinv = mean(cu(ct)./sigma).*rand(1);
    lb = [zeros(n,1);-Inf];
    ub = [ones(n,1);Inf];
    nlcon = [];
    theta_exp0 = [gamma0;cinv];

    [val,~,~] = fmincon(@(theta)obj(theta,beta,delta,r(ct,:),lambda_f_obs(ct),W_obs(ct,:)./sigma,J_obs(ct)./sigma, meanprices_fuel(ct,:),W(ct,:)),theta_exp0,[],[],[],[],lb,ub,nlcon,options);
    
    x_exp(ct,1) = val(1:n);
    x_exp(ct,2) = val(n+1);
    
    gamma_e = x_exp(ct,1);
    c_inv = x_exp(ct,2);
    
    t1 = (1-gamma_e)./(1-beta.*delta.*(1-gamma_e.*lambda_f_obs(ct)));
    t1 = t1.*(beta.*delta.*c_inv+(1-beta.*delta).*r(ct,:));

    t2 = gamma_e.*(1-beta.*delta.*(1-lambda_f_obs(ct)))./(1-beta.*delta.*(1-gamma_e.*lambda_f_obs(ct)));
    t2 = t2.*(J_obs(ct)-W_obs(ct,:))./sigma;
    meanprices_sim(ct,:) = t1+t2;
    
end

gamma_e = x_exp(:,1);
c_inv = x_exp(:,2);


%% Set up potential market 
%-------- First we compute the number of potential entrants. 

path_production = '';

ngp = 17;

W = ones(M,M);

%% Compute the exporters' value function 

d = freights(2:end,:)-(delta.*(freights(1:end-1,:)-matches(1:end-1,:)));   % new exporters observed in each market. Note that we can recover the
                                                                           % number of new exporters entering the market, d, from the steady state
                                                                           % equation f_(t+1) = d_t + delta*(f_t - m_t).




E = max(d(:,1:15),[],1)';

d = mean(d,1)';

G_unc_99 = zeros(M,M+1);

G_unc_99(:,1) = 1-(d./E);                                                  % potential exporters who decide not to enter

dd = repmat(d./E,[1 M]);

g = Gdest.*dd;                                                             % this describes the chosen destination 
                                                                           % of the potential exporters who decide to enter
g(eye(size(g))==1) = nan;

G_unc_99(:,2:end) = g;

meanprices_sim = zeros(M,M);            % given the estimates, we compute the simulated prices

for i = 1:M

    t1 = (1-gamma_e(i))./(1-beta.*delta.*(1-gamma_e(i).*lambda_f_obs(i)));
    t1 = t1.*(beta.*delta.*c_inv(i)+(1-beta.*delta).*r(i,:));

    t2 = gamma_e(i).*(1-beta.*delta.*(1-lambda_f_obs(i)))./(1-beta.*delta.*(1-gamma_e(i).*lambda_f_obs(i)));
    t2 = t2.*(J_obs(i)-W_obs(i,:))./sigma;
    meanprices_sim(i,:) = t1+t2;

end

Ue = -c_inv+lambda_f_obs.*(r- meanprices_sim); % and the corresponding exporters' value function
Ue = Ue./(1-beta.*delta.*(1-lambda_f_obs));

%% Inversion to recover the exporters entry costs
% We finally estimate the exporter entry costs, K, which capture both the cost of production, as well as
% any export costs beyond shipping prices. 
% To reduce noise, we cluster the entry costs into ngp groups.
% This is achieved by first estimating the entrire MxM matrix
% of entry costs, K_init. Then, we run a k-means algorithm to group
% the routes (ij) and re-estimate the entry costs separately for each 
% group

K_init = zeros(M,M);
    
for i=1:M

    ccp = log(G_unc_obs(i,2:end))-log(G_unc_obs(i,1));

    K_init(i,:) = -ccp + Ue(i,:);
end

ccp = zeros(M,M);
for i=1:M
    ss = min(G_unc_obs(i,:));
    ccp(i,:) = log(G_unc_obs(i,2:end)./ss)-log(G_unc_obs(i,1)./ss);
end


if ngp<sum(~isnan(K_init(:)))
    opts = statset('Display','off');
    [gp,~,sumd] = kmeans(K_init(:),ngp,'Replicates',100000,'Options',opts,'MaxIter',500000);
    gp(isnan(gp)) = ngp+1;

    X = zeros(length(ccp(:)),ngp);
    X = X(gp~=ngp+1,:);
    y = -ccp(:)+Ue(:);
    y = y(gp~=ngp+1);
    W = m_init.*Gdest_obs;
    W = ones(M,M);
    W = W(gp~=ngp+1);
    reg = fitlm(X,y,'Intercept',false,'Weights',W(:));
    K0 = reg.Coefficients.Estimate;
    
    K0 = [K0(1:ngp);0];
    K_obs = reshape(K0(gp),M,M);
    K_obs(eye(M)==1) = nan;
end
